#ifndef EM_PORT_API
#	if defined(__EMSCRIPTEN__)
#       include <emscripten/emscripten.h>
#		if defined(__cplusplus)
#			define EM_PORT_API(rettype) extern "C" rettype EMSCRIPTEN_KEEPALIVE
#		else
#			define EM_PORT_API(rettype) rettype EMSCRIPTEN_KEEPALIVE
#		endif
#	else
#		if defined(__cplusplus)
#			define EM_PORT_API(rettype) extern "C" rettype
#		else
#			define EM_PORT_API(rettype) rettype
#		endif
#	endif
#endif


#include <iostream>
#include <string>
#include <vector>
#include <map>
#include <tuple>
#include "math.h"

#include<detect_earthwork.h>


using namespace std;


extern  "C" {

class Point3f
{
public:
    float x;
    float y;
    float z;
    Point3f(float x0, float y0, float z0);
    ~Point3f();
};

Point3f::Point3f(float x0, float y0, float z0)
{
    x = x0;
    y = y0;
    z = z0;
}

Point3f::~Point3f()
{
}


int PtInPolygon(Point3f p, std::vector<Point3f>& ptPolygon)
{
	int nCross = 0;
	for (int i = 0; i < ptPolygon.size(); i++)
	{
		Point3f p1 = ptPolygon[i];
		Point3f p2 = ptPolygon[(i + 1) % ptPolygon.size()];

		if (p1.y == p2.y)
			continue;
		if (p.y < std::min(p1.y, p2.y))
			continue;
		if (p.y >= std::max(p1.y, p2.y))
			continue;

		float x = (p.y - p1.y) * (p2.x - p1.x) / (p2.y - p1.y) + p1.x;

		if (x > p.x)
		{
			nCross++;
		}

	}
	if (fabs(nCross % 2 - 1) < 0.01)
	{
		return 1;
	}
	else
	{
		return 0;
	}
}


EM_PORT_API(void) clip_points_with_volume(float * p_values, int p_value_count, float * poly_values, int poly_value_count, float * select_values, int * select_value_count, float * result)
{
	cout << "p_value_count" << p_value_count << "poly_value_count" << poly_value_count << endl;
	
	float max_dist = 100000000.0;
	float dist_eplison = 0.0001;
	//float grid_size = 0.3;
	float grid_size=0.015;
    float ref_elevation = -max_dist;

    int pts_count = p_value_count / 3;
    int poly_pts_count = poly_value_count / 3;
    std::vector<Point3f> poly;
	cout << "init poly" << endl;
    for (int i = 0; i < poly_pts_count;  i++)
    {
		cout << poly_values[3 * i] << ", " << poly_values[3 * i + 1] << ", " << poly_values[3 * i + 2] << endl;
        poly.push_back(Point3f(poly_values[3 * i], poly_values[3 * i + 1], 0));
        if (poly_values[3 * i + 2] > ref_elevation)
        {
          ref_elevation = poly_values[3 * i + 2];
        }
    }

	cout << "start selecting pts" << endl;
    select_value_count[0] = 0;
	int select_id = 0;
	// float min_x = max_dist, min_y = max_dist, max_x = -max_dist, max_y = -max_dist;
    float min_x =poly_values[0], min_y =poly_values[1], max_x = poly_values[0], max_y =poly_values[1];
	
    for (int i = 0; i < pts_count; i++)
	{
		if (i < 5)
		{
			cout << p_values[3 * i] << ", " << p_values[3 * i + 1] << ", " << p_values[3 * i + 2] << endl;
		}

		if (p_values[3 * i + 2] < ref_elevation - dist_eplison && PtInPolygon(Point3f(p_values[3 * i], p_values[3 * i + 1], 0), poly))
		{	
            select_values[3 * select_id] = p_values[3 * i]; // X
            select_values[3 * select_id + 1] = p_values[3 * i + 1]; // Y
            select_values[3 * select_id + 2] = p_values[3 * i + 2]; // Z

			if (p_values[3 * i] < min_x)
			{
				min_x = p_values[3 * i];
			}
			if (p_values[3 * i] > max_x)
			{
				max_x = p_values[3 * i];
			}
			if (p_values[3 * i + 1] < min_y)
			{
				min_y = p_values[3 * i + 1];
			}
			if (p_values[3 * i + 1] > max_y)
			{
				max_y = p_values[3 * i + 1];
			}
			
			select_id += 1;
			select_value_count[0] += 3;
		}
    }
	
	cout << "min_x: " << min_x << " max_x: " << max_x << " min_y: " << min_y << " max_y: " << max_y << endl;
	cout << "select_count: "  << select_id << endl;

	cout << "start init grid elevation: " << endl;
	// init grid elevation
	map<pair<int, int>, float> grid_elevation;
	int step_x = ceil((max_x - min_x) / grid_size);
	int step_y = ceil((max_y - min_y) / grid_size);
	cout << "step_x: " << step_x << ", step_y: " << step_y << endl;

	for (size_t i = 0; i < step_x; i++)
	{
		for (size_t j = 0; j < step_y; j++)
		{
			grid_elevation.insert(make_pair(make_pair(i, j), max_dist));
		}
	}

	cout << "get grid real elevation" << endl;
	// get grid real elevation
	for (size_t k = 0; k < select_id; k++)
	{
		int id_x = floor((select_values[3 * k] - min_x) / grid_size);
		int id_y = floor((select_values[3 * k + 1] - min_y) / grid_size);
		if (select_values[3 * k + 2] < grid_elevation[make_pair(id_x, id_y)])
		{
			grid_elevation[make_pair(id_x, id_y)] = select_values[3 * k + 2];
		}
	}

    int sum=0;
    for (size_t i = 0; i < step_x; i++)
    {
        for (size_t j = 0; j < step_y; j++)
        {
            if(grid_elevation[make_pair(i,j)]==max_dist)
            {
                sum++;
            }
            //grid_elevation.insert(make_pair(make_pair(i, j), max_dist));
        }
    }

    printf("once%d %lf\n",sum,sum/(double)(step_y*step_x));
	
	// calc volume
	cout << "calc volume" << endl;
	result[0] = 0.0;
	int valid_grid_count = 0;
	for (size_t i = 0; i < step_x; i++)
	{
		for (size_t j = 0; j < step_y; j++)
		{
			float real_elevation = grid_elevation[make_pair(i, j)];
			if (real_elevation < ref_elevation - dist_eplison)
			{
				float delta_height = ref_elevation - real_elevation;
				float cube_volume = grid_size * grid_size * delta_height;
				result[0] += cube_volume;
				valid_grid_count += 1;
			}
		}
	}

	cout << "valid_grid_count: " << valid_grid_count << endl;
	cout << "result: " << result[0] << endl;
}


EM_PORT_API(void) calc_earthwork_residual(float * bim_values, int bim_value_count, float * ew_values, int ew_value_count, float * poly_values, int poly_value_count,
                                                                                                float * selected_bim_values, int *selected_bim_value_count, float * selected_ew_values, int *selected_ew_value_count, float *residual_volume, float* bim_total_volume, float *residual_ratio)
{
    // log
    cout << bim_values[0] << ", " << bim_values[1] << ", " << bim_values[2] << endl;
    cout << ew_values[0] << ", " << ew_values[1] << ", " << ew_values[2] << endl;
    cout << "bim_value_count" << bim_value_count << "ew_value_count" << ew_value_count <<  "poly_value_count" << poly_value_count << endl;

    float max_dist = 100000000.0;
    float dist_eplison = 0.0001;
    float change_epsilon = 0.15;
    float earthwork_expand_ratio = 1.33;
    float grid_size = 0.3;
    float ref_elevation = -max_dist;

    int poly_pts_count = poly_value_count / 3;
    std::vector<Point3f> poly;
    cout << "init poly" << endl;
    float min_x = max_dist, min_y = max_dist, max_x = -max_dist, max_y = -max_dist;
    for (int i = 0; i < poly_pts_count;  i++)
    {
		cout << poly_values[3 * i] << ", " << poly_values[3 * i + 1] << ", " << poly_values[3 * i + 2] << endl;
        poly.push_back(Point3f(poly_values[3 * i], poly_values[3 * i + 1], 0));
        // update ref_elevation
        if (poly_values[3 * i + 2] > ref_elevation)
        {
            ref_elevation = poly_values[3 * i + 2];
        }

        // update min-max
        if (poly_values[3 * i] < min_x)
        {
            min_x = poly_values[3 * i];
        }
        if (poly_values[3 * i] > max_x)
        {
            max_x = poly_values[3 * i];
        }
        if (poly_values[3 * i + 1] < min_y)
        {
            min_y = poly_values[3 * i + 1];
        }
        if (poly_values[3 * i + 1] > max_y)
        {
            max_y = poly_values[3 * i + 1];
        }
    }

    cout << "ref_elevation: " << ref_elevation << endl;
    cout << "min_x: " << min_x << " max_x: " << max_x << " min_y: " << min_y << " max_y: " << max_y << endl;

    // init grid elevation
    cout << "start init grid elevation: " << endl;
    map<pair<int, int>, float> bim_elevation, earthwork_elevation;
    int step_x = ceil((max_x - min_x) / grid_size);
    int step_y = ceil((max_y - min_y) / grid_size);
    cout << "step_x: " << step_x << ", step_y: " << step_y << endl;

    for (size_t i = 0; i < step_x; i++)
    {
        for (size_t j = 0; j < step_y; j++)
        {
            bim_elevation.insert(make_pair(make_pair(i, j), max_dist));
        }
    }

    // copy bim_elevation to earthwork_elevation
	copy(bim_elevation.begin(), bim_elevation.end(), inserter(earthwork_elevation, earthwork_elevation.begin()));

    // update bim_elevation
    cout << "update bim_elevation" << endl;
    int bim_pt_count = bim_value_count / 3;
    float bim_max_elevation = -max_dist;
    for (int i = 0; i < bim_pt_count; i++)
    {
        if (PtInPolygon(Point3f(bim_values[3 * i], bim_values[3 * i + 1], 0), poly))
        {
            int id_x = floor((bim_values[3 * i] - min_x) / grid_size);
            int id_y = floor((bim_values[3 * i + 1] - min_y) / grid_size);

            if (bim_values[3 * i + 2] < bim_elevation[make_pair(id_x, id_y)])
            {
                bim_elevation[make_pair(id_x, id_y)] = bim_values[3 * i + 2];
            }

            // update bim max elevation
            if (bim_values[3 * i + 2] > bim_max_elevation)
            {
                bim_max_elevation = bim_values[3 * i + 2];
            }
        }
    }

    cout << "bim_max_elevation: "<< bim_max_elevation << endl;

    // update earthwork_elevation
    cout << "update earthwork_elevation" << endl;
    int ew_pt_count = ew_value_count / 3;
    for (int i = 0; i < ew_pt_count; i++)
    {
        if (PtInPolygon(Point3f(ew_values[3 * i], ew_values[3 * i + 1], 0), poly))
        {
            int id_x = floor((ew_values[3 * i] - min_x) / grid_size);
            int id_y = floor((ew_values[3 * i + 1] - min_y) / grid_size);

            if (ew_values[3 * i + 2] < earthwork_elevation[make_pair(id_x, id_y)])
            {
                earthwork_elevation[make_pair(id_x, id_y)] = ew_values[3 * i + 2];
            }
        }
    }
    
    // calc selected earthwork
    cout << " calc selected earthwork" << endl;
    int varied_grid_index = 0;
    residual_volume[0] = 0.0;
    float bim_volume = 0.0;
    for (size_t i = 0; i < step_x; i++)
    {
        for (size_t j = 0; j < step_y; j++)
        {
            float cur_bim_elevation = bim_elevation[make_pair(i, j)];
            float cur_earthwork_elevation = earthwork_elevation[make_pair(i, j)];
            float real_x = i * grid_size + min_x, real_y = j * grid_size + min_y;
            if (cur_bim_elevation < max_dist - dist_eplison && cur_earthwork_elevation < max_dist - dist_eplison && cur_earthwork_elevation >  cur_bim_elevation + change_epsilon)
            {
                // get selected bim cloud and selected earthwork cloud.
                selected_bim_values[3 * varied_grid_index] = real_x;
                selected_bim_values[3 * varied_grid_index + 1] = real_y;
                selected_bim_values[3 * varied_grid_index + 2] = cur_bim_elevation;
                selected_ew_values[3 * varied_grid_index] = real_x;
                selected_ew_values[3 * varied_grid_index + 1] = real_y;
                selected_ew_values[3 * varied_grid_index + 2] = cur_earthwork_elevation;

                float cube_volume = grid_size * grid_size * (cur_earthwork_elevation - cur_bim_elevation);
                // cube_volume *= earthwork_expand_ratio;
                residual_volume[0] += cube_volume;

                varied_grid_index += 1;
            }
            
            // calc bim volume
            if (cur_bim_elevation < max_dist - dist_eplison && bim_max_elevation > cur_bim_elevation + dist_eplison)
            {
                float cube_bim_volume = grid_size * grid_size * (bim_max_elevation - cur_bim_elevation);
                bim_volume += cube_bim_volume;
            }
        }
    }

    selected_bim_value_count[0] = 3 * varied_grid_index;
    selected_ew_value_count[0] = 3 * varied_grid_index;

    cout << "residual_volume: " << residual_volume[0] << endl;
    cout << "bim_volume: " << bim_volume << endl;

    // calc rest earthwork ratio
    residual_ratio[0] = residual_volume[0] / bim_volume;
    bim_total_volume[0] = bim_volume;

    bim_elevation.clear();
    earthwork_elevation.clear();

}


// cloud a is earlier than cloud b.
EM_PORT_API(void) detect_earthwork_change(float * values_A, int count_A, float * values_B, int count_B, float * poly_values, int poly_value_count,
                                                                  float * changed_values_A, int *changed_count_A, float * changed_values_B, int *changed_count_B, float *changed_volume)
{
    // log
    cout << values_A[0] << ", " << values_A[1] << ", " << values_A[2] << endl;
    cout << values_B[0] << ", " << values_B[1] << ", " << values_B[2] << endl;
    cout << "count_A" << count_A << "count_B" << count_B <<  "poly_value_count" << poly_value_count << endl;

    float max_dist = 100000000.0;
    float dist_eplison = 0.0001;
    float change_epsilon = 0.15;
    float earthwork_expand_ratio = 1.33;
    float grid_size = 0.3;
    float ref_elevation = -max_dist;

    int poly_pts_count = poly_value_count / 3;
    std::vector<Point3f> poly;
    cout << "init poly" << endl;
    float min_x = max_dist, min_y = max_dist, max_x = -max_dist, max_y = -max_dist;
    for (int i = 0; i < poly_pts_count;  i++)
    {
		cout << poly_values[3 * i] << ", " << poly_values[3 * i + 1] << ", " << poly_values[3 * i + 2] << endl;
        poly.push_back(Point3f(poly_values[3 * i], poly_values[3 * i + 1], 0));
        // update ref_elevation
        if (poly_values[3 * i + 2] > ref_elevation)
        {
            ref_elevation = poly_values[3 * i + 2];
        }

        // update min-max
        if (poly_values[3 * i] < min_x)
        {
            min_x = poly_values[3 * i];
        }
        if (poly_values[3 * i] > max_x)
        {
            max_x = poly_values[3 * i];
        }
        if (poly_values[3 * i + 1] < min_y)
        {
            min_y = poly_values[3 * i + 1];
        }
        if (poly_values[3 * i + 1] > max_y)
        {
            max_y = poly_values[3 * i + 1];
        }
    }

    cout << "ref_elevation: " << ref_elevation << endl;
    cout << "min_x: " << min_x << " max_x: " << max_x << " min_y: " << min_y << " max_y: " << max_y << endl;

    // init grid elevation
    cout << "start init grid elevation: " << endl;
    map<pair<int, int>, float> grid_elevation_A, grid_elevation_B;
    int step_x = ceil((max_x - min_x) / grid_size);
    int step_y = ceil((max_y - min_y) / grid_size);
    cout << "step_x: " << step_x << ", step_y: " << step_y << endl;

    for (size_t i = 0; i < step_x; i++)
    {
        for (size_t j = 0; j < step_y; j++)
        {
            grid_elevation_A.insert(make_pair(make_pair(i, j), max_dist));
        }
    }

    // copy grid_elevation_A to grid_elevation_B
	copy(grid_elevation_A.begin(), grid_elevation_A.end(), inserter(grid_elevation_B, grid_elevation_B.begin()));

    // update grid_elevation_A
    cout << "update grid_elevation_A" << endl;
    int pts_A_count = count_A / 3;
    for (int i = 0; i < pts_A_count; i++)
    {
        if (values_A[3 * i + 2] < ref_elevation - dist_eplison && PtInPolygon(Point3f(values_A[3 * i], values_A[3 * i + 1], 0), poly))
        {
            int id_x = floor((values_A[3 * i] - min_x) / grid_size);
            int id_y = floor((values_A[3 * i + 1] - min_y) / grid_size);

            if (values_A[3 * i + 2] < grid_elevation_A[make_pair(id_x, id_y)])
            {
                grid_elevation_A[make_pair(id_x, id_y)] = values_A[3 * i + 2];
            }
        }
    }

    // update grid_elevation_B
    cout << "update grid_elevation_B" << endl;
    int pts_B_count = count_B / 3;
    for (int i = 0; i < pts_B_count; i++)
    {
        if (values_B[3 * i + 2] < ref_elevation - dist_eplison && PtInPolygon(Point3f(values_B[3 * i], values_B[3 * i + 1], 0), poly))
        {
            int id_x = floor((values_B[3 * i] - min_x) / grid_size);
            int id_y = floor((values_B[3 * i + 1] - min_y) / grid_size);

            if (values_B[3 * i + 2] < grid_elevation_B[make_pair(id_x, id_y)])
            {
                grid_elevation_B[make_pair(id_x, id_y)] = values_B[3 * i + 2];
            }
        }
    }
    
    // calc earthwork changes
    cout << "calc earthwork changes" << endl;
    int varied_grid_index = 0;
    changed_volume[0] = 0.0;
    for (size_t i = 0; i < step_x; i++)
    {
        for (size_t j = 0; j < step_y; j++)
        {
            float elevation_A = grid_elevation_A[make_pair(i, j)];
            float elevation_B = grid_elevation_B[make_pair(i, j)];
            float real_x = i * grid_size + min_x, real_y = j * grid_size + min_y;
            if (elevation_A < ref_elevation - dist_eplison && elevation_B < ref_elevation - dist_eplison &&  fabs(elevation_A - elevation_B)  > change_epsilon)
            {
                // get changed cloud A and changed cloud B
                changed_values_A[3 * varied_grid_index] = real_x;
                changed_values_A[3 * varied_grid_index + 1] = real_y;
                changed_values_A[3 * varied_grid_index + 2] = elevation_A;
                changed_values_B[3 * varied_grid_index] = real_x;
                changed_values_B[3 * varied_grid_index + 1] = real_y;
                changed_values_B[3 * varied_grid_index + 2] = elevation_B;

                float cube_volume = grid_size * grid_size * (elevation_A - elevation_B);
                if (cube_volume > dist_eplison)
                {
                    cube_volume *= earthwork_expand_ratio;
                }
                changed_volume[0] += cube_volume;

                varied_grid_index += 1;
            }
        }
    }

    changed_count_A[0] = 3 * varied_grid_index;
    changed_count_B[0] = 3 * varied_grid_index;

    cout << "changed_volume: " << changed_volume[0] << endl;
    grid_elevation_A.clear();
    grid_elevation_B.clear();

}


// int main(){
//     // EM_ASM(console.log('你好，Emscripten！'));
//     std::cout<< "hello world!" <<std::endl;
//     return EXIT_SUCCESS;
// }


}